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ABSTRACT 


In realistic models of the crust and the upper mantle where vertical gradients in wave speeds and major interfaces 
are present, propagation of Pn, an important seismic phase at regional distances, involves complex effects of 
interference. Such effects can result in wave-trains whose frequency-contents and amplitudes vary with distance in 
counter-intuitive ways. In this new project, we are investigating the propagation of Pn beneath western China using 
data from Hi-CLIMB (An Integrated Study of the Himalayan-Tibetan Continental Lithosphere during Mountain 
Building). This experiment is one of the largest broadband seismic experiments to date, with more than 210 
deployments at close station-spacing of 3-8 km over a distance of 800 km. The linear array is complemented by a 
regional array of comparable aperture, producing an unprecedented dataset for Eurasia. For the current project, 
which Just began in May of 2008, we are organizing a dataset for regional seismic events recorded by the 
Hi-CLIMB arrays and will show examples of long seismic profiles over apertures of 500 km from several different 
azimuths. We are also investigating methods of modeling and inversion based on Gaussian beams (GB), which offer 
several distinct advantages. First, GB modeling can be applied both for interference waves, which have caustics, and 
for pure head waves. Second, GB modeling can handle laterally varying media, an important aspect that cannot be 
investigated by standard methods such as reflectivity. Third, GB is computationally efficient, suitable for analyzing 
large datasets. We will illustrate our approaches with data from an earlier experiment using explosions. Based on 
numerous events recorded by the Hi-CLIMB array, we will investigate frequency-dependent propagation of Pn over 
a large region in western China using GB and also finite-difference schemes. Attributes of seismic wave-trains, 
including arrival times, amplitudes of signal-envelopes, and instantaneous pulse frequencies will be used to 
constrain how structures in the crust and the upper mantle affect the propagation of Pn. The objective is to achieve 
self-consistent models of Pn propagation in western China which are free from assumptions such as a 
frequency-independent factor for geometric spreading. To this end, our results should advance efforts in isolating 
effects of frequency-dependent propagation from those of pure-inelastic attenuation (Q), leading to improved 
methodologies for discrimination and yield estimates at regional distances. 
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OBJECTIVES 

In this new project, which started in May 2008, we are analyzing the propagation of regional Pn waves recorded by 
Hi-CLIMB (An Integrated Study of the Himalayan-Tibetan Continental Lithosphere during Mountain Building) in 
Tibet. Considering the dense station-spacing (3 to 8 km) and large aperture (800 km) of the Hi-CLIMB seismic 
experiment, our objective is the use these densely recorded data to construct a self-consistent model of Pn 
propagation over a large region in western China. In turn, a better understanding of frequency-dependent effects of 
propagation from data inversion/modeling will improve frequency dependent attenuation models of regional Pn 
phases, leading to a clearer separation of inelastic from elastic effects in a laterally varying crust and upper mantle, 
eventually eliminating the assumption of a simplified, frequency-independent geometric spreading component. 

To account for the intricate nature of Pn propagation, we are applying the Gaussian beam (GB) method, an 
asymptotic approach. GB is stable at caustics and can be used for the diving waves, pure head-waves, and 
interference waves by choosing suitable beam parameters (Nowack and Stacy, 2002). A major advantage of GB is 
that it works well in laterally varying media. In contrast, alternatives such as the reflectivity method, are restricted to 
one-dimensional media. When compared with purely numerical simulations for laterally varying media, such as the 
finite difference method, GB is computationally efficient in facilitating modeling and inversion of large 
datasets—tasks in which a large number of forward calculations are required. However, for small scale scattering, 
hybrid methods or full finite difference modeling may still be required. As a practical measure, essential features of 
Pn waveforms will be distilled into seismic attributes such as arrival times, envelopes of wave amplitudes, and 
instantaneous pulse frequencies for modeling and inversion. 

RESEARCH ACCOMPLISHED 

The Hi-CLIMB seismic array was deployed between 2002 and 2005, in three phases, over 210 broadband sites at a 
station-spacing of 3-8 km (for the main, linear array). The main, north-south trending linear array is complemented 
by an east-west expanding regional array. The overall the aperture is about 800 km by 800 km (Figure 1). A unique 
feature of Hi-CLIMB is that it is the only seismic profile that crosses the highest portions of the Himalayas, linking 
baseline measurements of the Indian shield, in a pristine state prior to its involvement in active deformation, with 
data in the interior of Tibet. There are abundant, high quality data from numerous events at regional distances from 
the array. Figure 2 shows a preliminary plot of regional events with magnitudes greater than 4.6 over the time period 
of the Hi-CLIMB deployment. 

Figure 3 shows a preliminary seismic profile from a single event with a magnitude of 4.6 located to the west of the 
array (located by the star in Figure 2) and recorded by central and northern portions of the Hi-CLIMB array. Clear 
arrivals of both Pn and Pg phases are observed to epicentral distances of more than 550 km, with Pn overtaking Pg 
as the first arrival at distances near 300 km. The first arrivals are marked by the dotted line, and the Pn arrival 
becomes a first arrival after about 325 km. By comparing the signals marked by box A near 350 km with a dominant 
frequency of about .5 Hz with box B near 560 km where the pulse frequency is higher at about .7 Hz, this increase in 
dominant frequency with distance suggests complex interference effects for Pn propagation in this region. A 
preliminary profile from a second event with a magnitude of 4.8 with a location just south of the middle and 
northern portions of the array is shown in Figure 4. Here, the Pn becomes the first arrival after about 290 km. The 
character of the Pn is more complicated for this profile than for the profile for the first event, and this suggests that 
there is significant lateral variability in the region that will need to be addressed in this study. 

To interpret first arrival Pn data recorded from regional events by the Hi-CLIMB array, the Gaussian beam method 
will be used for larger scale features in the model. The Gaussian beam method (GB) is an asymptotic approach for 
the computation of high-frequency wavefields in laterally varying media by the summation of paraxial beams, each 
surrounding a central ray-path and having a Gaussian distribution in amplitude away from the central ray (Popov, 
1982; Cerveny et al., 1982; Nowack and Aki, 1984). Reviews of the method are given by Babich and Popov (1989) 
and more recently by Nowack (2003) and Popov (2002). There are several advantages in using GB. First, individual 
beam components have no singularities along their paths. Thus the summation of beams is also regular everywhere, 
in contrast to ray methods. Second, GB naturally introduces some smoothing and therefore is not as sensitive to 
model parameterizations as ray-methods. Third, GB does not require the added computation of two-point ray 
tracing. 
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New approaches are being developed to perform two parameter beam summations based on initial position and 
angle of individual beam components. The advantage of these new approaches is that a wider range of beam 
parameters can be used for the each individual beam solution, but the trade-off is in speed of computation 
(Lugara et al., 2003; Nowack, 2003; Nowack et al., 2006; Nowack et al., 2007a; Nowack, 2008). 

(b) 




Figure 1. Map showing geologic setting and configuration of Project Hi-CLIMB. (a) Grey-scale topographic 
map showing overall setting of the Himalayan-Tibetan collision zone. Two open arrows show 
ground velocity with respect to stable Eurasia. The box shows the region depicted in panel (b). Thin, 
solid curves show major geological boundaries, (b) Topographic map showing the full configuration 
of the Hi-CLIMB broadband seismic array. Notice a large aperture of 800 km by 800 km and the 
dense spacing of stations along the main, linear array. In addition to the IRIS/PASSCAL pool, 
instruments also came from several other participating institutions (see legend). 



m 


Figure 2. A color-coded topographic map, showing epicenters of large to moderate-sized earthquakes 

(mt 4.6 or above; open red circles) recorded by the Hi-CLIMB array (purple triangles). The asterisk 
marks the epicenter of a Tibetan earthquake (/Wft~4.8) which serves as the seismic source of the 
preliminary profile shown in Figure 4. 


To efficiently describe essential features of a complex waveform, it is useful to identify key attributes of the data. 
Procedures for extracting such attributes include complex trace analysis, and the general approach of wavelet 
analysis (e.g., Daubechies, 1992; Mallat, 1999). The approach described here is the former, in which an analytic 
signal is constructed from a seismogram. From the analytic signal, its envelope and instantaneous frequency can be 
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determined. Applications of complex trace analysis to seismic data have been described by Taner et al. (1979), and 
an application to estimate seismic attenuation by matching observed values of instantaneous frequency was 



Distcxe (tTi) 


Figure 3. An example of seismic profile from a single event (marked by the star in Figure 2) recorded by 

central and northern portions of the Hi-CLIMB main array, showing clear arrivals of Pn (marked 
by the dotted curve) which become the first arrival at distances beyond about 325 km. Notice the 
dominant frequency of Pn actually can increase with distance: Compare signals marked by box A, 
about 0.5 Hz near 350 km; and box B, about 0.7 Hz near 550 km. Solid circles mark first arrivals at 
selected seismograms. The time-scale is reduced by a wave speed of 8.0 km/s. 



Oistcnce (krn) 

Figure 4. A second example of a seismic profile from a single event with a magnitude of 4.8 located directly 
south of the middle and northern sections of the Hi-CLIMB array. The Pn becomes the first arrival 
at distances beyond about 290 km. The time-scale is reduced by a wave speed of 8.0 km/s. 

performed by Matheney and Nowack (1995). It is important to recall that the instantaneous frequency is generally 
not equivalent to the spectral frequency. However, an appropriate weighted average of the instantaneous frequency 
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converges to an average of positive spectral frequencies. Figure 5 shows an example of results from this procedure. 
For a seismic phase of interest, one can extract amplitudes from its envelope, timing of the seismic phase, and 
instantaneous frequencies as a function of time. For estimating seismic attenuation, one can match instantaneous 
frequencies between observed and attenuated reference pulses (Matheney and Nowack, 1995). More general 
inversions of all seismic attributes for selected seismic phases are also feasible, including previous work by Nowack 
and Matheney (1997) and Matheney et al. (1997). 


OBSERVED TRACE 



WEIGHTED & DAMPED INSTANTANEOUS FREQUENCY 



T me (lec) 

Figure 5. An example of extracting attributes from seismic data. A) Comparison of observed seismogram 

(solid trace) and its envelope calculated based on the Hilbert transform of data (dashed curve). The 
peak amplitude of the Hilbert envelope for the first pulse is marked. B) Damped and smoothed 
instantaneous frequency as a function of time (from Matheney and Nowack, 1995). 


Pn is not a simple head-wave but is comprised of complex sets of interfering waves when there is a positive velocity 
gradient with respect to depth below the Moho interface (Cerveny and Ravindra, 1971; Hill, 1971; Menke and 
Richards, 1980). These interfering waves can simply result from the spherical nature of the earth (Hill, 1973; Sereno 
and Given, 1990; Yang et al., 2007), an effect equivalent to a velocity gradient in the corresponding flattened model. 
How observed amplitudes of refracted arrivals change over distance supports the interpretation that these 


Vclocit> (iradtcni Mixicl 
V docii> (im t) 



Pure r>i\ing Wa\e C ^ 




Figure 6. Illustration showing how a positive gradient in wave velocity beneath the Moho can affect the 
propagation of Pn, A) The model for P-wave speeds. The linear gradient beneath the Moho can 
simply arise from the spherical nature of the earth. B) Paths of the primary, diving rays (pure head 
waves are not shown). C) Paths of the first set of sub-Moho interference wave. D) Paths of the 
second set of sub-Moho interference wave (from Nowack and Stacy, 2002). 
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arrivals are comprised of diving and interference waves, instead of only a pure head wave (Braile and Smith, 1975). 
Figure 6 shows an example of how a positive upper mantle gradient can effect the propagation of Pn waves by 
forming a set of interference waves at the Moho (from Nowack and Stacy, 2002). 

In the last several years, there has been a debate on the effects of small scale scattering in the upper mantle on the 
propagation of teleseismic Pn. Ryberg et al. (1995, 2000) and Tittgemeyer et al. (1996) favored a model in which 
scattering in the upper mantle can account for many of the characteristics of teleseismic Pn. In contrast, Morozov et 
al. (1998) and Morozov and Smithson (2000) advocated a model with a positive upper mantle gradient forming an 
interference head wave along with crustal scattering to explain teleseismic Pn. However, the recent analysis of 
Nielson and Thybo (2003) supports the model that the Pn phase travels as an upper mantle whispering gallery and 
that the origin of its long coda at teleseismic distances is due to crustal scattering. For our study, we will model the 
first arrival Pn phase as an interference wave at regional distances using Gaussian beam modeling. At greater 
distances, a hybrid approach or the direct use of finite difference modeling may be required to model the 
characteristics of teleseismic Pn as a combination of upper mantle gradients and small scale scattering in the crust or 
upper mantle. 

Nowack and Stacy (2002) compared expansions of interference waves using Gaussian beams and the reflectivity 
method, which is restricted to one-dimensional velocity structures. The fact that GB expansions are non-singular at 
caustics is essential for multiple underside reflections ,which are associated with caustics and pseudo-caustics, and 
for these cases finite-sized beams are used. Meanwhile, broad beams are used to represent wide-angle reflections, 
including the pure head-wave (Nowack and Aki, 1984). In both cases, beams are specified by effective planar 
phase-fronts at the receiver, providing stable results during summations. 

Figure 7 shows comparisons of results from GB and reflectivity, where the reflectivity modeling includes multiples 
from the free surface which were not included in the GB modeling. In general, these two methods yield comparable 
results: Notice the good agreement between the two methods in key crustal and Pn seismic attributes at all distances. 
In Fig. 7C, the envelope amplitudes of the first arriving pulses are shown following the procedure of Mathaney and 
Nowack (1995). The gap near 125 km corresponds to the cross-over distance between the direct crustal arrival and 
Pn. Fig. 7D shows the instantaneous pulse frequencies for the first arrival. At distances less than 100 km, 
instantaneous pulse frequencies are for the direct, crustal arrivals. At greater distances, the values are for Pn. For this 
particular velocity model, instantaneous frequencies (fi) initially drop just beyond the cross-over distance, mainly 
reflecting properties of the pure head-wave. Subsequently,starts recovering as diving and interference waves 
begin to emerge.// continues to increase with distance and eventually overshoots beyond values of direct arrivals 
seen at close-in distances, an effect resulting from interference of different sets of underside reflections below the 
Moho. This result is confirmed by analyzing centroid-frequencies of the amplitude spectra of the relevant pulses 
(Nowack and Stacy, 2002). A similar trend in/ was reported by Cerveny and Ravindra (1971) using asymptotic 
analysis. These pure-elastic effects in regional body-waves must be accounted for in order to obtain meaningful 
estimates of attenuation. 

It should be emphasized again that unlike reflectivity, GB is valid for laterally varying media. This is an essential 
feature to study regional body-waves in areas such as western China where laterally varying structure is most likely 
to be prevalent in the lithosphere. To this end, we will further validate GB in laterally varying earth models by 
comparing with other methods which generates complete wavefields, such as the finite difference method 
(Virieux, 1988). The advantage of using the Gaussian beam method is primarily in its computational speed, a 
necessary feature for inversion or modeling of massive datasets where a large number of forward simulations must 
be performed. 

The discussed methodologies to be used for Hi-CLIMB were previously applied to data from the Deep 
Probe/SAREX project by Stacy and Nowack (2002). SAREX (Southern Alberta Refraction Experiment) was 
designed to provide detailed information on crustal structures along the Canadian portion of the combined Deep 
Probe/SAREX profile (Clowes et al., 2002; Gorman et al., 2002). To interpret the observations, Stacy and Nowack 
(2002) computed synthetic seismograms using GB; whose attributes were then compared with those extracted from 
the data. The resulting features in the amplitudes and the instantaneous frequencies were found to be more 
complicated than that shown in Figure 7. In addition, a model for Q was simultaneously determined, resulting in a 
self-consistent representation of both elastic and inelastic effects of wave propagation over a large range of 
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distances. When modeling Pn propagation beneath Hi-CLIMB in western China, it is anticipated that more 
complicated Pn propagation characteristics will also be found. 
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Figure 7 Synthetic seismograms and associated seismic attributes from the simple layer over a gradient 

model. A) Synthetic seismograms using the reflectivity method. B) Synthetic seismograms using the 
Gaussian beam method. C) Amplitude of the envelope of the first arriving pulse as a function of 
distance. Notice that there is a sudden drop in amplitudes near distances of 100 km. Results from 
Gaussian beam modeling (pluses) and those from reflectivity (squares) are similar. D) Instantaneous 
frequency of the first arrival pulse as a function of distance. Beyond distance of about 120 km, 
instantaneous frequency actually increases with increasing distance. Again, results from Gaussian 
beam modeling (pluses) and those from reflectivity (squares) are similar. The source pulse is a Gabor 
wavelet with a dominant frequency of 2.77 Hz. To isolate elastic effects of wave propagation, a 
constant Q value of 2000 was assumed (from Nowack and Stacy, 2002). 

A major initial task of the project is to first make a comprehensive database of regional events recorded by the 
Hi-CLlMB array. The database will include full waveforms, instrument responses, a catalog of events and other 
important information such as uncertainties of locations and focal mechanisms. Meissner et al. (2004) used 
earthquake data recorded by the INDEPTH 111 array in central Tibet to study regional phases and concluded that 
errors in epicenters of up to 15 km are negligible for Pn and would only slightly affect Pg at close-in distances. So 
errors in epicenter do not appear to be a major concern. Nevertheless, we plan to quickly assess such errors by 
comparing event catalogs based on global (Engdahl, et al., 1998) and regional datasets, both of which are available 
in the public domain (the latter is available at the China Earthquake Networks Center, CENC). In terms of focal 
depths, it has been known for some time that most events are shallow in Tibet, between 5 to 10 km (Molnar and 
Chen, 1983; Chen and Molnar, 1983; Langin et al., 2003). Nonetheless, if necessary, relocations will be performed 
using a combination of regional data from the Hi-CLIMB array and teleseismic data. One possibility is to adopt a 
simplified relocation approach, as used by Zhu et al. (2006) for regional events in Tibet. 

If necessary, for larger events, nih of 5.5 and above, we can obtain precise results of source parameters from the 
inversion of waveforms recorded at teleseismic distances. Similarly, for smaller events, their source parameters can 
be retrieved from inversion of regional waveforms (Zhu et al., 2006). Both procedures are routinely performed by 
Chen’s group (e.g., Kao et al., 2002; Chen and Yang, 2004). Meanwhile, there are a number of models reported in 
the literature for the crust and upper mantle of Tibet (e.g., Hearn et al., 2004; Liang et al., 2004, 2006; Sun and 
Toksoz, 2006, Li et al., 2006), sufficient for constructing preliminary background models. Also, a recent study by 
Phillips et al. (2007) used Pn travel-times to invert for upper mantle velocities and gradients in the region. 
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Once a suitable set of events, with precise locations of epicenters, depths, and focal mechanisms, is assembled, we 
will bring together a starting velocity model from published results, including those from recent Pn travel-time 
tomography and seismic refraction studies. In addition, we will incorporate our own results from migration imaging 
using teleseismic data from the Hi-CLIMB array (Nowack et al., 2007b). We shall then build a large number of 
seismic profiles based on these chosen events. In the first example (Figure 3), an increase of frequency for Pn with 
propagation distance is already apparent (cf Boxes A and B in Figure 3), suggesting that the Pn from this region is 
being influenced by interference effects. The next step is to extract seismic attributes from each seismic profile so 
that quantitative analysis and modeling can be carried out. 

The modeling of these attributes involves matching observed envelope-amplitudes and instantaneous pulse 
frequencies with those calculated, mainly using Gaussian beams. To minimize effects of seismic source-time 
functions, one strategy is to match pulses observed at close-in distances with theoretical ones and normalize the rest 
of the data accordingly. Alternatively, the normalization can be performed by deconvolution of all data with an 
appropriate reference signal. The reference can be a P-pulse observed at teleseismic distance, if the event is large 
enough (w/, above ~5 or 5.5). Otherwise, a pulse observed at close-in distances could be appropriate. Once the 
seismic attributes are normalized, we will then perform systematic modeling using Gaussian beams. The final 
product will be several laterally varying velocity and attenuation models over a large region of western China, since 
numerous regional earthquake from different azimuths were well-recorded by the Hi-CLIMB array (Figure 2). 

CONCLUSIONS AND RECOMMENDATIONS 


For the first year of this new project, we are building a comprehensive database of regional events recorded by the 
Hi-CLIMB array. The starting point of this database will be relevant global and regional catalogs of seismicity, 
followed by relocation of selected events and inversion for source depths and focal mechanisms, if required. The 
major deliverable will be the database, including the catalog of events with precise locations of epicenters, focal 
depths, and focal mechanisms. Further validation of the Gaussian beam method is being conducted to efficiently and 
precisely calculate Pn propagation in laterally varying crustal and upper mantle structures. We will carefully 
compare results from different methods, including those from finite difference algorithms. The comparisons will 
also include extraction of relevant seismic attributes. 

In the second phase of the project, we will construct seismic profiles across the Hi-CLIMB array for a selected set of 
relocated regional events. Normalization and deconvolution of seismic profiles by various reference signals in order 
to remove effects of the source-time function will be performed. Extraction of key seismic attributes will then be 
performed from observed seismic profiles, and detailed modeling/inversion of seismic attributes using Gaussian 
beams will be conducted in order to build self-consistent models for propagation of Pn beneath the Hi-CLIMB array 
in western China. 
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